Set Up

library(tidyverse)
library(here)
library(phyloseq)
library(vegan)
theme_set(theme_bw())
filt_fun <- function (x, min_reads = 100, min_samples = 5) {
    (sum(x) > min_reads) & (sum(x > 0) > min_samples)
}
ps <- readRDS(here('data','following_study','ps_rarefied.rds')) 
sample_data(ps)[is.na(sample_data(ps))] <- 'NA'
# transform data into proportion
ps.prop <- ps %>% transform_sample_counts(function(x) x/sum(x))
# double check proportion transform
sample_sums(ps.prop) %>% summary()
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1
sam <- data.frame(sample_data(ps))
max.core <- parallel::detectCores()

Alpha Diversity

alpha.diversity <- data.frame(sample_data(ps),estimate_richness(ps))

Household Effect

ggplot(alpha.diversity,aes(x = as.numeric(Household), y = Shannon, group = Household)) +
    geom_point() + geom_line() + xlab('Household')

anova(lm(Shannon~Household, data = alpha.diversity))
## Analysis of Variance Table
## 
## Response: Shannon
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Household 48 8.9301 0.18604  1.5916 0.05425 .
## Residuals 49 5.7277 0.11689                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Epileptic Effect

ggplot(alpha.diversity,aes(x = Epileptic.or.Control, y = Shannon)) +
    geom_boxplot() + geom_jitter(height = 0, width = 0.25)

Paired t test

epileptic <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
control <- alpha.diversity %>% filter(Epileptic.or.Control == 'Control')
t.test(epileptic$Shannon, control$Shannon, paired = TRUE)
## 
##  Paired t-test
## 
## data:  epileptic$Shannon and control$Shannon
## t = -0.22089, df = 48, p-value = 0.8261
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.1556564  0.1248409
## sample estimates:
## mean difference 
##     -0.01540773

ANOVA with Household effect added

anova(lm(Shannon~Household + Epileptic.or.Control, data = alpha.diversity))
## Analysis of Variance Table
## 
## Response: Shannon
##                      Df Sum Sq  Mean Sq F value  Pr(>F)  
## Household            48 8.9301 0.186045  1.5607 0.06329 .
## Epileptic.or.Control  1 0.0058 0.005816  0.0488 0.82612  
## Residuals            48 5.7219 0.119206                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Breed Effect

ggplot(alpha.diversity) +
    geom_point(aes(x = Breed.Group..1., y = Shannon, colour = Breed.Group..1.)) +
    facet_wrap(~Epileptic.or.Control) + 
    theme(axis.text.x = element_blank(), axis.ticks.x.bottom = element_blank())

# here breed group with na is removed 
anova(lm(Shannon~Household + Breed.Group..1., data = alpha.diversity %>% filter(Breed.Group..1. != 'NA')))
## Analysis of Variance Table
## 
## Response: Shannon
##                 Df Sum Sq Mean Sq F value  Pr(>F)  
## Household       45 7.4046 0.16455  1.5646 0.07819 .
## Breed.Group..1.  4 1.2267 0.30668  2.9160 0.03339 *
## Residuals       39 4.1016 0.10517                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Drug Effect

alpha.diversity.epi <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
ggplot(alpha.diversity.epi, aes(x = Pheno.Y.N, y = Shannon)) +
    geom_boxplot() + geom_jitter(height = 0, width = 0.25)

t.test(Shannon~Pheno.Y.N, data = alpha.diversity.epi)
## 
##  Welch Two Sample t-test
## 
## data:  Shannon by Pheno.Y.N
## t = 0.98407, df = 33.649, p-value = 0.3321
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##  -0.1316909  0.3787797
## sample estimates:
##  mean in group No mean in group Yes 
##          2.856864          2.733320

Sex Effect

ggplot(alpha.diversity, aes(x = Sex, y = Shannon)) +
    geom_boxplot() + geom_jitter(height = 0, width = 0.25)

anova(lm(Shannon~Household + Sex, data = alpha.diversity))
## Analysis of Variance Table
## 
## Response: Shannon
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Household 48 8.9301 0.18604  1.6309 0.04673 *
## Sex        1 0.2522 0.25218  2.2107 0.14360  
## Residuals 48 5.4755 0.11407                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(Shannon~Sex, data = alpha.diversity)
## 
##  Welch Two Sample t-test
## 
## data:  Shannon by Sex
## t = -0.83555, df = 88.411, p-value = 0.4057
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -0.2223635  0.0907199
## sample estimates:
## mean in group F mean in group M 
##        2.792318        2.858140

Age Effect

ggplot(alpha.diversity, aes(x = as.numeric(Age..months.)/12, y = Shannon)) +
    geom_point() + xlab('Age')
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

lm(Shannon~as.numeric(Age..months.), data = alpha.diversity) %>% summary()
## Warning in eval(predvars, data, env): NAs introduced by coercion
## 
## Call:
## lm(formula = Shannon ~ as.numeric(Age..months.), data = alpha.diversity)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90685 -0.30999  0.01342  0.28351  0.75221 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              2.745594   0.091949  29.860   <2e-16 ***
## as.numeric(Age..months.) 0.000838   0.001114   0.752    0.454    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3838 on 94 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.005986,   Adjusted R-squared:  -0.004588 
## F-statistic: 0.5661 on 1 and 94 DF,  p-value: 0.4537

Beta Diversity

Test for Effect

plot_ord <- function(data, factor, method, distance){
    data.ord <- ordinate(data, method = method, distance = distance)
    p <- plot_ordination(data, data.ord, color = factor)
    p <- p + stat_ellipse(type = "t",geom = "polygon",alpha = 0)
    p <- p + ggtitle(str_c(factor,method,distance, sep = ' - '))
    print(p)
}
permanova <- function(data, formula, method, permutations=1e4, strata = NULL, core = max.core){
    dist.matrix <- phyloseq::distance(data, method=method)
    df <- data.frame(sample_data(data))
    model <- as.formula(paste0('dist.matrix~', formula))
    if (!is_null(strata)) {strata <- df[,strata]}
    result <- adonis2(model,
                      data = df,
                      permutations=permutations,
                      strata = strata,
                      parallel = core)
    return(result)
}
permdisp <- function(data, group, method, permutations=1e4, pairwise = FALSE, core = max.core){
    dist.matrix <- phyloseq::distance(data, method=method)
    df <- data.frame(sample_data(data))
    beta.disp <- betadisper(dist.matrix, group = df[,group])
    result <- permutest(beta.disp, permutations = permutations, pairwise = pairwise, type = 'centroid')
    return(result)
}

Household Effect

permanova(ps.prop, 'Household', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F    Pr(>F)    
## Household 48  12.2165 0.68881 2.2596 9.999e-05 ***
## Residual  49   5.5192 0.31119                     
## Total     97  17.7358 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Epileptic Effect

Ordination

plot_ord(ps.prop, 'Epileptic.or.Control','MDS','bray')

plot_ord(ps.prop, 'Epileptic.or.Control','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2235 -- in the phylogenetic tree in the data you provided.

plot_ord(ps.prop, 'Epileptic.or.Control','NMDS','bray')
## Run 0 stress 0.2086595 
## Run 1 stress 0.2041723 
## ... New best solution
## ... Procrustes: rmse 0.05570532  max resid 0.4115003 
## Run 2 stress 0.2049336 
## Run 3 stress 0.2224345 
## Run 4 stress 0.2036331 
## ... New best solution
## ... Procrustes: rmse 0.03306726  max resid 0.305915 
## Run 5 stress 0.204184 
## Run 6 stress 0.2031278 
## ... New best solution
## ... Procrustes: rmse 0.01942513  max resid 0.0845663 
## Run 7 stress 0.2098507 
## Run 8 stress 0.2099703 
## Run 9 stress 0.2217262 
## Run 10 stress 0.2258699 
## Run 11 stress 0.2040622 
## Run 12 stress 0.2229979 
## Run 13 stress 0.2110726 
## Run 14 stress 0.2042057 
## Run 15 stress 0.2264617 
## Run 16 stress 0.229214 
## Run 17 stress 0.2194904 
## Run 18 stress 0.2031184 
## ... New best solution
## ... Procrustes: rmse 0.01466582  max resid 0.08599131 
## Run 19 stress 0.2053133 
## Run 20 stress 0.2155131 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      8: no. of iterations >= maxit
##     12: stress ratio > sratmax

plot_ord(ps.prop, 'Epileptic.or.Control','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV75 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1827303 
## Run 1 stress 0.1948716 
## Run 2 stress 0.1826111 
## ... New best solution
## ... Procrustes: rmse 0.02702614  max resid 0.2334849 
## Run 3 stress 0.1901819 
## Run 4 stress 0.1827303 
## ... Procrustes: rmse 0.027013  max resid 0.2336186 
## Run 5 stress 0.1878658 
## Run 6 stress 0.1882704 
## Run 7 stress 0.1910674 
## Run 8 stress 0.1821596 
## ... New best solution
## ... Procrustes: rmse 0.03709371  max resid 0.2358942 
## Run 9 stress 0.1813226 
## ... New best solution
## ... Procrustes: rmse 0.04120045  max resid 0.2334379 
## Run 10 stress 0.1904447 
## Run 11 stress 0.1870262 
## Run 12 stress 0.1879809 
## Run 13 stress 0.1869531 
## Run 14 stress 0.1809398 
## ... New best solution
## ... Procrustes: rmse 0.04006268  max resid 0.2347613 
## Run 15 stress 0.1896303 
## Run 16 stress 0.1883852 
## Run 17 stress 0.1896765 
## Run 18 stress 0.1966149 
## Run 19 stress 0.189808 
## Run 20 stress 0.1831356 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

PERMANOVA Test

permanova(ps.prop, 'Epileptic.or.Control', 'bray', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                      Df SumOfSqs      R2      F Pr(>F)
## Epileptic.or.Control  1   0.1554 0.00876 0.8487 0.1039
## Residual             96  17.5803 0.99124              
## Total                97  17.7358 1.00000
permanova(ps.prop, 'Epileptic.or.Control', 'wunifrac', strata = 'Household')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2981 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                      Df SumOfSqs      R2      F  Pr(>F)  
## Epileptic.or.Control  1  0.02437 0.01471 1.4336 0.05139 .
## Residual             96  1.63226 0.98529                 
## Total                97  1.65664 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMDISP

permdisp(ps.prop, 'Epileptic.or.Control', 'bray')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.00167 0.0016733 0.1117  10000 0.7396
## Residuals 96 1.43851 0.0149844
permdisp(ps.prop, 'Epileptic.or.Control', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1119 -- in the phylogenetic tree in the data you provided.
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.0002514 0.00025138 2.0869  10000 0.1522
## Residuals 96 0.0115640 0.00012046

Breed Effect

a table of # of breed. frequency table

(data.frame(sample_data(ps.prop))$Breed.Group..1.) %>% table()
## .
##          Herder              NA Pointer Spaniel       Retriever     Scent Hound 
##              20               9              26              26               5 
##     Sight Hound        Sled Dog         Terrier 
##               3               6               3
plot_ord(ps.prop, 'Breed.Group..1.','MDS','bray')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plot_ord(ps.prop, 'Breed.Group..1.','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2583 -- in the phylogenetic tree in the data you provided.
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plot_ord(ps.prop, 'Breed.Group..1.','NMDS','bray')
## Run 0 stress 0.2086595 
## Run 1 stress 0.2055186 
## ... New best solution
## ... Procrustes: rmse 0.05947923  max resid 0.4053172 
## Run 2 stress 0.2217881 
## Run 3 stress 0.2086136 
## Run 4 stress 0.2223882 
## Run 5 stress 0.2042025 
## ... New best solution
## ... Procrustes: rmse 0.01751769  max resid 0.1099105 
## Run 6 stress 0.2041345 
## ... New best solution
## ... Procrustes: rmse 0.01516255  max resid 0.08075781 
## Run 7 stress 0.2111371 
## Run 8 stress 0.218612 
## Run 9 stress 0.2030952 
## ... New best solution
## ... Procrustes: rmse 0.03875405  max resid 0.3146495 
## Run 10 stress 0.2121138 
## Run 11 stress 0.2043533 
## Run 12 stress 0.2041224 
## Run 13 stress 0.2134545 
## Run 14 stress 0.2041461 
## Run 15 stress 0.2260212 
## Run 16 stress 0.2087215 
## Run 17 stress 0.2249867 
## Run 18 stress 0.2053203 
## Run 19 stress 0.2050026 
## Run 20 stress 0.2046693 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      3: no. of iterations >= maxit
##     17: stress ratio > sratmax
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plot_ord(ps.prop, 'Breed.Group..1.','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV731 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1910295 
## Run 1 stress 0.1947317 
## Run 2 stress 0.1939297 
## Run 3 stress 0.1892023 
## ... New best solution
## ... Procrustes: rmse 0.04324271  max resid 0.2330391 
## Run 4 stress 0.1915156 
## Run 5 stress 0.1943251 
## Run 6 stress 0.1884852 
## ... New best solution
## ... Procrustes: rmse 0.03014957  max resid 0.2137321 
## Run 7 stress 0.1919418 
## Run 8 stress 0.1943292 
## Run 9 stress 0.1877478 
## ... New best solution
## ... Procrustes: rmse 0.06971062  max resid 0.2875466 
## Run 10 stress 0.1928138 
## Run 11 stress 0.1900625 
## Run 12 stress 0.1932147 
## Run 13 stress 0.1910529 
## Run 14 stress 0.195628 
## Run 15 stress 0.1929301 
## Run 16 stress 0.1967131 
## Run 17 stress 0.1891158 
## Run 18 stress 0.1968316 
## Run 19 stress 0.1999649 
## Run 20 stress 0.1943933 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     18: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

PERMANOVA

# same breed group v.s. diff breed group in household
# check if data are paired within household
if (!identical(sample_data(ps)$Household[seq(1,98,2)],
               sample_data(ps)$Household[seq(1,98,2)+1])) stop('data is not paired')

same.breed.group <- table(sample_data(ps)$Household, sample_data(ps)$Breed.Group..1.) %>% 
    apply(1, function(x) any(x == 2)) # if dog's grouping are same with in each house
same.breed.group;sum(same.breed.group)
##     1    10    11    12    13    14    15    16    17    18    19     2    20 
## FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
##    21    22    23    24    25    26    27    28    29     3    30    31    32 
##  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE 
##    33    34    35    36    38    39     4    40    41    42    43    44    45 
##  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
##    47    48    49     5    50    51     6     7     8     9 
##  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
## [1] 39

Here we see 39 out of 49 household have dogs in same breed group

permanova(ps.prop, 'Household + Breed.Group..1.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F    Pr(>F)    
## Household       48  12.2165 0.68881 2.3397 9.999e-05 ***
## Breed.Group..1.  6   0.8417 0.04746 1.2896   0.07309 .  
## Residual        43   4.6775 0.26374                     
## Total           97  17.7358 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1666 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F    Pr(>F)    
## Household       48  0.82865 0.67272 2.2283 9.999e-05 ***
## Breed.Group..1.  6  0.07000 0.05683 1.5059   0.05389 .  
## Residual        43  0.33314 0.27045                     
## Total           97  1.23179 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMANOVA Test by Epileptic and Control Group

Here we

ps.control <- ps.prop %>% subset_samples(Epileptic.or.Control == 'Control')
ps.epileptic <- ps.prop %>% subset_samples(Epileptic.or.Control == 'Epileptic')
perm.breed.bray.ctl <- permanova(ps.control, 'Breed.Group..1.', 'bray')
perm.breed.wunif.ctl <- permanova(ps.control, 'Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2158 -- in the phylogenetic tree in the data you provided.
perm.breed.bray.ctl;perm.breed.wunif.ctl
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs    R2      F Pr(>F)  
## Breed.Group..1.  7   1.7534 0.194 1.4098 0.0295 *
## Residual        41   7.2851 0.806                
## Total           48   9.0385 1.000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df  SumOfSqs      R2      F Pr(>F)  
## Breed.Group..1.  7 0.0033166 0.21306 1.5858 0.0386 *
## Residual        41 0.0122498 0.78694                
## Total           48 0.0155664 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perm.breed.bray.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'bray')
perm.breed.wunif.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2069 -- in the phylogenetic tree in the data you provided.
perm.breed.bray.epi;perm.breed.wunif.epi
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F Pr(>F)
## Breed.Group..1.  7   1.3363 0.15645 1.0863   0.26
## Residual        41   7.2055 0.84355              
## Total           48   8.5418 1.00000
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F Pr(>F)
## Breed.Group..1.  7  0.21498 0.17043 1.2033  0.178
## Residual        41  1.04644 0.82957              
## Total           48  1.26142 1.00000

Combined p value

\[X^2=-2 \sum_{i=1}^k \ln \left(p_i\right)\]

bray.p <- c(perm.breed.bray.ctl$`Pr(>F)`[1], perm.breed.bray.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(bray.p)), df = 4, lower.tail = FALSE)
## [1] 0.04501879
wunif.p <- c(perm.breed.wunif.ctl$`Pr(>F)`[1], perm.breed.wunif.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(wunif.p)), df = 4, lower.tail = FALSE)
## [1] 0.0410838

With using the weighted-unifrac distance, there’s at least one of the null hypothesis is rejected.

Drug Effect

Here we are only focusing on epileptic dogs.

plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','bray')

plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2547 -- in the phylogenetic tree in the data you provided.

plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','bray')
## Run 0 stress 0.2020935 
## Run 1 stress 0.1992304 
## ... New best solution
## ... Procrustes: rmse 0.04145127  max resid 0.1509698 
## Run 2 stress 0.2182018 
## Run 3 stress 0.1994009 
## ... Procrustes: rmse 0.06905381  max resid 0.3691135 
## Run 4 stress 0.2034522 
## Run 5 stress 0.2002274 
## Run 6 stress 0.2117299 
## Run 7 stress 0.2094952 
## Run 8 stress 0.2029557 
## Run 9 stress 0.2082424 
## Run 10 stress 0.2087795 
## Run 11 stress 0.2081331 
## Run 12 stress 0.2041039 
## Run 13 stress 0.2025866 
## Run 14 stress 0.2106487 
## Run 15 stress 0.2086674 
## Run 16 stress 0.2145508 
## Run 17 stress 0.2136685 
## Run 18 stress 0.2061055 
## Run 19 stress 0.2084818 
## Run 20 stress 0.2100444 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      3: no. of iterations >= maxit
##     17: stress ratio > sratmax

plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV3027 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1803735 
## Run 1 stress 0.18135 
## Run 2 stress 0.1940171 
## Run 3 stress 0.1831955 
## Run 4 stress 0.192244 
## Run 5 stress 0.181868 
## Run 6 stress 0.1814495 
## Run 7 stress 0.1856995 
## Run 8 stress 0.1827014 
## Run 9 stress 0.1859298 
## Run 10 stress 0.1803533 
## ... New best solution
## ... Procrustes: rmse 0.097442  max resid 0.3798744 
## Run 11 stress 0.1785945 
## ... New best solution
## ... Procrustes: rmse 0.1067785  max resid 0.5529312 
## Run 12 stress 0.1873117 
## Run 13 stress 0.1856336 
## Run 14 stress 0.186715 
## Run 15 stress 0.1780467 
## ... New best solution
## ... Procrustes: rmse 0.09321676  max resid 0.2595204 
## Run 16 stress 0.1857093 
## Run 17 stress 0.1849503 
## Run 18 stress 0.1945162 
## Run 19 stress 0.1873817 
## Run 20 stress 0.1782206 
## ... Procrustes: rmse 0.0993426  max resid 0.270433 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

PERMANOVA

permanova(ps.epileptic, 'Pheno.Y.N','bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F Pr(>F)
## Pheno.Y.N  1   0.2211 0.02588 1.2488  0.183
## Residual  47   8.3207 0.97412              
## Total     48   8.5418 1.00000
permanova(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV794 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F Pr(>F)
## Pheno.Y.N  1  0.01552 0.01978 0.9482 0.4564
## Residual  47  0.76925 0.98022              
## Total     48  0.78477 1.00000

Dispersion

permdisp(ps.epileptic, 'Pheno.Y.N','bray')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.02029 0.020293 1.5629  10000 0.2151
## Residuals 47 0.61028 0.012985
permdisp(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2646 -- in the phylogenetic tree in the data you provided.
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)  
## Groups     1 0.0004601 0.00046008 4.9251  10000 0.0292 *
## Residuals 47 0.0043905 0.00009341                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# within epileptic dogs
table(sample_data(ps)$Seizure.Freedom..Y.N., sample_data(ps)$Epileptic.or.Control)
##      
##       Control Epileptic
##   NA       49         3
##   No        0        32
##   Yes       0        14
table(sample_data(ps)$Seizure.Control..Satisfactory.Unsatisfactory., sample_data(ps)$Epileptic.or.Control)
##     
##      Control Epileptic
##   NA      49         3
##   S        0        33
##   US       0        13

Sex Effect

Does health condition related to dog’s sex? Here we conduct a Chi-squared test.

# chisq test on epi male dog v.s. female.
# not within health condition
tb <- table(sample_data(ps)$Sex, sample_data(ps)$Epileptic.or.Control)
tb
##    
##     Control Epileptic
##   F      33        25
##   M      16        24
chisq.test(tb)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tb
## X-squared = 2.0698, df = 1, p-value = 0.1502

Ordination

plot_ord(ps, 'Sex','MDS','bray')

plot_ord(ps, 'Sex','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1364 -- in the phylogenetic tree in the data you provided.

plot_ord(ps, 'Sex','NMDS','bray')
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.242098 
## Run 1 stress 0.2453813 
## Run 2 stress 0.2402009 
## ... New best solution
## ... Procrustes: rmse 0.06907647  max resid 0.2808975 
## Run 3 stress 0.2497949 
## Run 4 stress 0.2425178 
## Run 5 stress 0.2450577 
## Run 6 stress 0.2467007 
## Run 7 stress 0.2470747 
## Run 8 stress 0.2524667 
## Run 9 stress 0.2505825 
## Run 10 stress 0.2536592 
## Run 11 stress 0.2438436 
## Run 12 stress 0.2468012 
## Run 13 stress 0.2552667 
## Run 14 stress 0.2526795 
## Run 15 stress 0.2630919 
## Run 16 stress 0.2425955 
## Run 17 stress 0.2442017 
## Run 18 stress 0.2600297 
## Run 19 stress 0.2606059 
## Run 20 stress 0.244812 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

plot_ord(ps, 'Sex','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2524 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1904303 
## Run 1 stress 0.1969604 
## Run 2 stress 0.1890081 
## ... New best solution
## ... Procrustes: rmse 0.03639876  max resid 0.1910498 
## Run 3 stress 0.18442 
## ... New best solution
## ... Procrustes: rmse 0.06855117  max resid 0.3274062 
## Run 4 stress 0.1843571 
## ... New best solution
## ... Procrustes: rmse 0.01999591  max resid 0.1783847 
## Run 5 stress 0.1911162 
## Run 6 stress 0.1856643 
## Run 7 stress 0.1856803 
## Run 8 stress 0.1980304 
## Run 9 stress 0.1902333 
## Run 10 stress 0.1901071 
## Run 11 stress 0.1886344 
## Run 12 stress 0.19008 
## Run 13 stress 0.1863481 
## Run 14 stress 0.1856805 
## Run 15 stress 0.1860512 
## Run 16 stress 0.1904701 
## Run 17 stress 0.1898851 
## Run 18 stress 0.1925002 
## Run 19 stress 0.1863481 
## Run 20 stress 0.1907652 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

PERMANOVA

permanova(ps.prop, 'Household + Sex', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F    Pr(>F)    
## Household 48  12.2165 0.68881 2.2373 9.999e-05 ***
## Sex        1   0.0589 0.00332 0.5174    0.9774    
## Residual  48   5.4604 0.30787                     
## Total     97  17.7358 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Sex', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1628 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F    Pr(>F)    
## Household 48  1.21161 0.67613 2.0949 9.999e-05 ***
## Sex        1  0.00201 0.00112 0.1669    0.9978    
## Residual  48  0.57836 0.32275                     
## Total     97  1.79198 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMDISP

permdisp(ps.prop, 'Sex', 'bray')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.01829 0.018289 1.2389  10000 0.2721
## Residuals 96 1.41725 0.014763
permdisp(ps.prop, 'Sex', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1999 -- in the phylogenetic tree in the data you provided.
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.000131 0.00013143 0.0899  10000 0.7675
## Residuals 96 0.140320 0.00146166

Age Effect

sam$Age..months. <- as.numeric(sam$Age..months.)
ggplot(sam) +
    geom_line(aes(x = as.numeric(Household), y = Age..months./12, group = Household)) + 
    geom_point(aes(x = as.numeric(Household), y = Age..months./12, colour = Epileptic.or.Control)) +
    xlab('Household') + ylab('Age in Year')

sam.epi <- sam %>% filter(Epileptic.or.Control == 'Epileptic')
sam.ctl <- sam %>% filter(Epileptic.or.Control == 'Control')
age.diff <- data.frame(Household = sam.epi$Household, Age.Diff = sam.epi$Age..months. - sam.ctl$Age..months.)
ggplot(age.diff) +
    geom_bar(aes(x = as.numeric(Household), y = Age.Diff/12), stat = 'identity') + 
    xlab('Household') + ylab('Age Difference in Year')

PERMANOVA

permanova(ps.prop, 'Household + Age..months.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##              Df SumOfSqs      R2      F    Pr(>F)    
## Household    48  12.2165 0.68881 2.3084 9.999e-05 ***
## Age..months. 19   2.2116 0.12470 1.0557     0.315    
## Residual     30   3.3076 0.18650                     
## Total        97  17.7358 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Age..months.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2330 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##              Df SumOfSqs      R2      F    Pr(>F)    
## Household    48 0.054452 0.67232 2.1861 9.999e-05 ***
## Age..months. 19 0.010971 0.13546 1.1127    0.2822    
## Residual     30 0.015567 0.19221                     
## Total        97 0.080990 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1